   #include "stdafx.h"
   
   using namespace std;
   
  Mod_PDF::Mod_PDF(Mod_in *iinp):Mod_out_handling(iinp){;}
  Mod_PDF::~Mod_PDF(){;}
int   Mod_PDF::compute(time_intervall &ti){}
  int Mod_PDF_Hydraulics_KA5::flag_soil_hydraulic_tabs_KA4;
  int Mod_PDF_Hydraulics_KA5::flag_soil_hydraulic_tabs_KA4_humus_zusatz;
  Par_Maps Mod_PDF_Hydraulics_KA5::soil_hydraulic_tabs_KA4;
  Par_Maps Mod_PDF_Hydraulics_KA5::soil_hydraulic_tabs_KA4_humus_zusatz;
  
  
   
   Mod_PDF_Hydraulics_KA5::Mod_PDF_Hydraulics_KA5(
		  map <string,int> &iIO_options,
		  Mod_in *iinp,
		  Mod_Soilphysics *imod_soilphysics,
		  Mod_wc *imod_wc):Mod_PDF(iinp)
	{
		 //1. switches setzen
		IO_options=iIO_options;
	
                
		//3. Paras_setzen
                flag_Mod_PDF_Hydraulics_already_set=0;
                if (flag_soil_hydraulic_tabs_KA4!=1)
                {
                    soil_hydraulic_tabs_KA4.read_Par_Maps("par","soil_hydraulic_tabs_KA4","par");
                    flag_soil_hydraulic_tabs_KA4=1;
                }
                 if (flag_soil_hydraulic_tabs_KA4_humus_zusatz!=1)
                {
                    soil_hydraulic_tabs_KA4.read_Par_Maps("par","soil_hydraulic_tabs_KA4_humus_zusatz","par");
                    flag_soil_hydraulic_tabs_KA4_humus_zusatz=1;
                }
		//4. Modelle verkn�pfen
		
		mod_soilphysics= imod_soilphysics;
		mod_wc=imod_wc;
		//5.outputs setzen
		 
		 define_field_of_states("Porosity",inp->horizons->size());
         define_field_of_states("FK_pf4.2",inp->horizons->size());
         define_field_of_states("FK_pf1.8",inp->horizons->size());
         define_field_of_states("FK_pf2.8pf",inp->horizons->size());
         define_field_of_states("FK",inp->horizons->size());
         define_field_of_states("KS",inp->horizons->size());
         define_field_of_states("William_b1",inp->horizons->size());
		 define_field_of_states("William_a1",inp->horizons->size());
	}
 Mod_PDF_Hydraulics_KA5::~Mod_PDF_Hydraulics_KA5(){;}  
 int   Mod_PDF_Hydraulics_KA5::compute(time_intervall &ti)
 {

	int y=0,z=0, decision,sand,silt,clay;
        string bd;
        int layer;
        timestamp start=ti.get_start();
        timestamp end=ti.get_end();
        for (layer=0;layer<inp->horizons->size();layer++)
        {
            if (flag_Mod_PDF_Hydraulics_already_set==1)
            {
                adjust_FC_for_groundwater_table(layer,mod_wc->get_states((string)"Grounwater_table_cm",start),start);
            }
            else
            {
            y=0,z=0;
            do
            {
                    //searching for the soiltype, determined by sand,silt,clay content
                    sand = (mod_soilphysics->get_states("sand",layer,start)>=soil_hydraulic_tabs_KA4.get_val("sand_low",y))
                            +(mod_soilphysics->get_states("sand",layer,start)<=soil_hydraulic_tabs_KA4.get_val("sand_up",y));
                    silt = (mod_soilphysics->get_states("silt",layer,start)>=soil_hydraulic_tabs_KA4.get_val("silt_low",y))
                            +(mod_soilphysics->get_states("silt",layer,start)<=soil_hydraulic_tabs_KA4.get_val("silt_up",y));
                    clay = (mod_soilphysics->get_states("clay",layer,start)>=soil_hydraulic_tabs_KA4.get_val("clay_low",y))
                            +(mod_soilphysics->get_states("clay",layer,start)<=soil_hydraulic_tabs_KA4.get_val("clay_up",y));
                    sand=sand+silt+clay;
                    y=y+1;

            }
            while (sand<6);
            y=y-1;
            if (mod_soilphysics->get_states("bd",layer,start)<= 1.45)
            {
                    bd="bd_low";//low bulkdensity
            }
            else
            {
                    if (mod_soilphysics->get_states("bd",layer,start)>=1.65)
                    {
                            bd="bd_med";//medium bulkdensity
                    }
                    else
                    {
                            bd="bd_high";// high bulkdensity
                    }
            }
            do
            {
                    decision=0;
                    decision=(soil_hydraulic_tabs_KA4.get_val("soiltyp_id",y)==soil_hydraulic_tabs_KA4_humus_zusatz.get_val("soiltyp_id",z));
                    decision=decision+(soil_hydraulic_tabs_KA4_humus_zusatz.get_val("humus_low",z)<=mod_soilphysics->get_states("SOC",layer,start)*1.72);
                    decision=decision+(soil_hydraulic_tabs_KA4_humus_zusatz.get_val("humus_high",z)>=mod_soilphysics->get_states("SOC",layer,start)*1.72);
                    if (decision==3)
                    {

                            this->set_states("Porosity",layer,end,
                                    ( soil_hydraulic_tabs_KA4.get_val("FK_"+bd,y)
                                    + soil_hydraulic_tabs_KA4.get_val("LK_"+bd,y)
                                    +soil_hydraulic_tabs_KA4_humus_zusatz.get_val("GPV",z) )/100);//wcs
                            this->set_states("FC_pf1.8",layer,end,
                                    ( soil_hydraulic_tabs_KA4.get_val("FK_"+bd,y)
                                    +soil_hydraulic_tabs_KA4_humus_zusatz.get_val("FK",z) )/100);//fc 1.8 pf
                            this->set_states("FC_pf2.5",layer,end,
                                    ( soil_hydraulic_tabs_KA4.get_val("FK_300_"+bd,y)
                                    +soil_hydraulic_tabs_KA4_humus_zusatz.get_val("FK",z) )/100);//fc fc pf2.5 starke Vereinfachung
                            this->set_states("FC_pf4.2",layer,end,
                                    ( soil_hydraulic_tabs_KA4.get_val("FK_"+bd,y)
                                    -soil_hydraulic_tabs_KA4.get_val("NFK_"+bd,y)
                                    +soil_hydraulic_tabs_KA4_humus_zusatz.get_val("GPV",z) 
                                    -soil_hydraulic_tabs_KA4_humus_zusatz.get_val("NFK",z)
                                    -+soil_hydraulic_tabs_KA4_humus_zusatz.get_val("LK",z) )/100);//fc fc pf4.2
                             this->set_states("FC_pf2.5",layer,end,
                                    ( soil_hydraulic_tabs_KA4.get_val("KF_LD_"+bd,y) )/100);//conductivity
                            //this->set_states("FC_pf1.8",ti.get_end(),soil [y] [12+bd]+humus[z][5])/100;// fc pf1.8
                            //this->set_states("FC_pf2.5",ti.get_end(),soil [y] [18+bd]+humus[z][5])/100;// fc pf2.5 starke Vereinfachung
                            //this->set_states("FC_pf4.2",ti.get_end(),soil [y] [12+bd]-soil[y][9+bd]+humus[z][6]-humus[z][3]-humus[z][4])/100;// fc pf4.2
                            //this->set_states("KS",ti.get_end(),soil [y] [15+bd]);// conductivity

                            if (mod_wc->get_states("Grounwater_table_cm",start)<100){
                                adjust_FC_for_groundwater_table(layer,mod_wc->get_states("Grounwater_table_cm",start),start);
                            }
                            else
                            {
                                    this->set_states("FC",end,(this->get_states("FC_pf1.8",end)+this->get_states("FC_pf2.5",end))/2);
                            }
                            //hydro_props [x][6][0]=(hydro_props [x][1][0]+hydro_props [x][2][0])/2;//Sch�tzwert f�r fieldcapacity
                    }
                    else
                    {
                            z=z+1;
                    }
                    //printf("%i  wcs %f  fc %f\n",z,humus [z][6],humus[z][5]);
            }
            while (decision<3);
          }
            compute_William_parameter(layer,end);
           
        }
	   return 0;
 }
 
 int Mod_PDF_Hydraulics_KA5::compute_William_parameter (int layer,timestamp &ts)
 {
     double theta_s=this->get_states("Porosity",layer,ts);
     double theta_r=this->get_states("FC_pf4.2",layer,ts);
     double h_s=5;
     double h_r=16000;
     
     double b1=(log(h_s)-log(h_r))/(log(theta_s)-log(theta_r));
     double a1=log(h_s)-b1*log(theta_s);
     
     this->set_states("William_b1",layer,ts,b1);
     this->set_states("William_a1",layer,ts,a1);
     return 0;
 }

 
 int Mod_PDF_Hydraulics_KA5::adjust_FC_for_groundwater_table(int layer,double gw_table,timestamp &ts)
 {
     
     // differenze of mid layer and gw_table
     double mid_layer;
     double diff;
     double FC;
     if(layer==0) mid_layer=inp->horizons->at(layer)/2;
     else mid_layer  =inp->horizons->at(layer-1)+(inp->horizons->at(layer)-inp->horizons->at(layer-1))/2;
     diff=gw_table-mid_layer;
     if (diff<=0)	// layer is below gw table
     {
         this->set_states("FC",layer,ts,this->get_states("Porosity",layer,ts));
     }
     else
     {
         //kleiner als 1.8
         if (diff<pow(1.8,10))
         {
             FC=0.5*(this->get_states("Porosity",layer,ts)+this->get_states("FC_pf1.8",layer,ts));
             this->set_states("FC",layer,ts,FC);
         }
          if ((diff<pow(2.27,10))&&(diff>=pow(1.8,10)))
         {
             FC=0.5*(this->get_states("FC_pf2.5",layer,ts)+this->get_states("FC_pf1.8",layer,ts));
             FC=0.5*(FC+this->get_states("FC_pf1.8",layer,ts));
             this->set_states("FC",layer,ts,FC);
         }
     }
     return 0;
 }